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The importance of the bulk viscosity of QCD in ultrarelativistic heavy-ion collisions 
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We investigate the consequences of a nonzero bulk viscosity coefficient on the transverse momen¬ 
tum spectra, azimuthal momentum anisotropy, and multiplicity of charged hadrons produced in 
heavy ion collisions at LHC energies. The agreement between a realistic 3D hybrid simulation and 
the experimentally measured data considerably improves with the addition of a bulk viscosity coeffi¬ 
cient for strongly interacting matter. This paves the way for an eventual quantitative determination 
of several QCD transport coefficients from the experimental heavy ion and hadron-nucleus collision 
programs. 


1. Introduction. Ultrarelativistic heavy ion collisions 
realized at the Relativistic Heavy Ion Collider (RHIC) 
and the Large Hadron Collider (LHC) are able to reach 
energies high enough to create and study the quark-gluon 
plasma (QGP), a novel state of nuclear matter where the 
quark and gluon degrees of freedom become manifest [T] . 
This hot and dense nuclear medium was found to behave 
like an almost perfect fluid, with one of the smallest shear 
viscosity to entropy density ratios, 77 /s, in nature [ 2 HS]. 
Currently, one of the main theoretical challenges in nu¬ 
clear physics is to model such collisions and extract from 
experiment the transport properties of this new phase of 
nuclear matter. 

Fluid-dynamical models have been highly successful in 
describing the production of hadrons in heavy ion colli¬ 
sions. The azimuthal momentum anisotropy of hadrons 
in particular has been shown to be a sensitive probe of the 
shear viscosity of the QGP [7] , and has been used repeat¬ 
edly to estimate this transport coefficient [S]. One limita¬ 
tion of this extraction procedure is the uncertainty asso¬ 
ciated with the early time dynamics of the collisions: the 
azimuthal momentum distribution of hadrons is known 
to be closely related to the initial shape of the medium 
[SHII]. Therefore, an accurate determination of the shear 
viscosity and other transport properties of QCD matter 
demands further improvements in the modeling of the 
earliest stages of the collisions. 

Recent improvements in modeling the early time dy¬ 
namics of heavy ion collisions [nma using the IP-Sat 
model of the nucleon wavefunction [M] followed by a clas¬ 
sical Yang-Mills evolution of the gluon fields [15] led to 
unprecedented success m in describing charged hadron 
azimuthal momentum distributions as characterized by 
their harmonic coefficients Vn {n = 2, 3,4,- • •). Further 
support for this initial state model, known as IP-Glasma, 
was provided by the remarkable agreement with data of 
its prediction for the event-by-event distributions of 
measured by the ATLAS collaboration m 

The same approach, however, had less success in de¬ 
scribing the full transverse momentum distribution of 
hadrons, showing clear tension with data in the low trans¬ 
verse momentum region m- In this letter we show that 


the inclusion of bulk viscosity, which was neglected in 
previous studies, can relieve this tension. In principle, 
the bulk viscosity of QCD matter should not be zero for 
the temperatures achieved at the RHIC and the LHC and 
it may become large enough to affect the evolution of the 
medium. In fact, simulations of heavy ion collisions that 
include the effect of bulk viscosity have already been per¬ 
formed [T5H25] and demonstrated that bulk viscosity can 
have a non-negligible effect on heavy ion observables. 

In addition to the early time description of the colli¬ 
sion provided by the IP-Glasma model, our calculations 
include a phase of hadronic re-scatterings after the hy¬ 
drodynamic evolution, implemented using the ultra rela¬ 
tivistic quantum molecular dynamics simulation UrQMD 
ESI [30]. Moreover, the intermediate fluid-dynamical evo¬ 
lution is resolved using a more complete version |31] of 
Israel-Stewart theory |32| that takes into account all the 
second order terms that couple the shear-stress tensor 
and bulk viscous pressure. This hybrid approach with 
IP-Glasma initial conditions is found to be capable of 
describing simultaneously the multiplicity and average 
transverse momentum of pions, kaons, and protons when 
a finite bulk viscosity, of the order C/s ~ 0.3, is included 
near the QCD phase transition region. Such a finite bulk 
viscosity also considerably reduces, by almost 50%, the 
value of the shear viscosity needed to describe the har¬ 
monic flow coefficients. 

2. Model. The initial state of the medium is deter¬ 
mined using the IP-Glasma model with the thermaliza- 
tion time set to tq = 0.4 fm. The system then evolves 
following the conservation law = 0, where the 

stress-energy tensor is composed of the ideal part 
= eu^u^ — A^'^Po{e) and the dissipative part = 
Here e is the local energy density, Po{£) is 
the thermodynamic pressure according to the equation of 
state, the fluid velocity, H the bulk viscous pressure, 
and the shear-stress tensor. We further introduced 
the projection operator — u^u'" onto the 3-D 

space orthogonal to the fluid velocity. The equation of 
state, Po{£), is the chemical equilibrium one taken from 
Ref. |33| . It is a parametrization of a lattice QCD calcu¬ 
lation matched onto a hadron resonance gas calculation 
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at lower temperatures. We assume that the baryon num¬ 
ber density and diffusion are zero at all space-time points 
and our metric convention is = diag(l, —1 — 1 — 1 ). 

The time-evolution equations satisfied by If and 
are relaxation-type equations derived from kinetic theory 
[311 [3S]. These are solved numerically within the MUSIC 
hydrodynamics simulation |36H38j . Explicitly, we solve 

+ =—C^O — , ( 1 ) 

+ ( 2 ) 

The above equations include all the nonlinear terms that 
couple bulk viscous pressure and shear-stress tensor and 
have recently been shown to be in good agreement with 
solutions of the 0-1-1 Anderson-Witting equation in the 
massive limit [39] and of the 1-1-1 Anderson-Witting equa¬ 
tion in the massless limit [101 El] • For the sake of simplic¬ 
ity, the transport coefficients rn, ^nn, Ahtt, Utt, V, 6-k-k, ‘Pt, 
Ttt-ki and A^n-n are fixed using formulas derived from the 
Boltzmann equation near the conformal limit [35]. The 
shear viscosity coefhcient is assumed to be proportional 
to the entropy density, i.e., 77 oc s. The bulk viscosity co¬ 
efficient employed is the same one introduced in Ref. [21] , 
which corresponds to a parametrization of calculations 
from Ref. [12] for the QGP phase and Ref. [H] for the 
hadronic phase. These two calculations are matched at 
Tc = 180 MeV and the value of C,/s at this tempera¬ 
ture is (^/s{Tc) « 0.3. This parametrization is plotted 
in Fig. 1 as the blue solid curve. The results shown in 
this letter have a small sensitivity on the value of (/s 
near the matching temperature, which can be doubled 
without leading to major modifications in our results. 



FIG. 1: (Color online) The bulk viscosity over entropy density 
parametrization used in our simulations as a function of T/Tc. 

At an isothermal hypersurface specified by the switch¬ 
ing temperature Tswitch, the simulation switches from 
a fluid-dynamical description to a transport description 


[44] ■ modeled using the UrQMD simulation. The mo¬ 
mentum distribution of hadrons at each hypersurface el¬ 
ement is calculated via the usual Cooper-Frye formalism 
pS] . The multiplicity of each hadron species is sampled 
assuming that every fluid element is a grand-canonical 
ensemble while the momentum of each hadron is obtained 
by sampling the momentum distribution using the rejec¬ 
tion method. We note that the Cooper-Frye formalism 
requires as an input the nonequilibrium momentum dis¬ 
tribution of each hadron inside the fluid elements. For 
the correction related to bulk viscous pressure, we employ 
the distribution derived from the Boltzmann equation us¬ 
ing the relaxation time approximation, as described in 
Ref. |46| . For the shear-stress tensor nonequilibrium cor¬ 
rection, we employ the usual ansatz obtained from the 
14-moment approximation [231EZ]- The details of how 
UrQMD is matched to MUSIC will be presented in an 
upcoming paper. 

We emphasize that the nonequilibrium corrections to 
the momentum distribution of hadrons at the moment 
of switching are still not completely understood from a 
theoretical point of view and represent a source of un¬ 
certainty in simulations of heavy ion collisions. However, 
the differential observables carry most of these uncertain¬ 
ties since they are more sensitive to the details of how 
the momentum of hadrons is distributed when convert¬ 
ing from a hydrodynamic to a transport description. For 
this reason, we fix all the free parameters of our model 
using pT“integrated observables. 

3. Results and Discussion. In our simulations, the 
value of the shear viscosity coefhcient is adjusted to pro¬ 
vide a good agreement with the integrated how harmonic 
coefficients, u„, up to n = 4. For the simulations that in¬ 
clude both bulk and shear viscosity, this procedure led 
to the value gjs = 0.095. For the simulations which in¬ 
clude only the shear viscosity, our baseline calculation 
is carried out with 77 /s = 0.16. The larger value of 77 /s 
compensates the reduction of momentum anisotropy due 
to the effect of the bulk viscosity. 

The pion and kaon multiplicity, and, to a lesser 
extent, their average px, are only mildly sensitive to 
the choice of switching temperature between the hydro- 
dynamic and UrQMD phases. Proton observables, on 
the other hand, do depend significantly on the choice of 
’^switch- The switching temperature used in the following 
calculations is fixed such that a good description of the 
proton multiplicity and average px is achieved for the 
simulation with both shear and bulk viscosities. This 
value is Tswitch = 145 MeV. 

In Figs. 2 (a), (b), and (c), we show the multiplic¬ 
ity, average transverse momentum of pions, kaons, and 
protons, and the integrated flow harmonics of charged 
hadrons, as a function of centrality class. The u„{2} co¬ 
efficients are calculated following the cumulant method 

[45] using the same px cuts employed by the ALICE col¬ 
laboration |49| . The multiplicity and average transverse 
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FIG. 2: (Color online) Multiplicity (a), average transverse momentum (b), and flow harmonic coefficients (c) as a function 
of centrality. The bands around the dashed lines show the effect of Tswitch on the observables. The points correspond to 
measurements by the ALICE collaboration |49l I50| . with bars denoting the experimental uncertainty. 


momentum are calculated without a lower px cut m- 
All resonances and hadrons included in UrQMD are con¬ 
sidered in our analyses and we neglect all weak decays. 
The solid curves correspond to the simulations that in¬ 
clude bulk and shear viscosities, while the dashed lines 
correspond to the calculations with only the shear vis¬ 
cosity. The band around the dashed curves shows how 
the results are modified when Tgwitch is varied from 135 
MeV to 165 MeV. For (px) and v„, the upper section of 
the band corresponds to the calculations with the lowest 
^switch while for multiplicity it corresponds to ones with 
the highest Tswitch- The points correspond to measure¬ 
ments by the ALICE collaboration [49l [50]. 

As expected, the simulations without bulk viscosity are 
still able to well describe the centrality dependence of 
the flow harmonic coefficients U2,3,4(2}. However, these 
calculations overestimate the (px) of pions, kaons, and 
protons by almost 30%. This happens because the IP- 
Glasma model gives rise to an initial state with large 
gradients of pressure and the subsequent fluid-dynamic 
expansion accordingly produces a significant radial flow. 
Therefore, in order to describe the data the transverse 
momentum of produced particles must be considerably 
reduced. 

Including hadronic re-scatterings by itself does not re¬ 
duce the (pt)) modifying mostly the intermediate px re¬ 
gion of the pion spectra Moreover, we can see 

from the bands around the dashed lines in Fig. 2 that 
increasing the switching temperature will not help fixing 
the multiplicity of pions, and is not enough to reproduce 
the correct values of (px)- Finally, reducing rj/s alone 
not only is unable to sufficiently suppress the (px), but 
also ends up destroying the good description of the flow 
harmonic coefficients. 

Including bulk viscosity leads to a suppression of (px) 
and can improve our description of the data. This is be¬ 


cause the bulk viscous pressure acts as a resistance to 
the expansion or compression of the fluid. In heavy ion 
collisions, the expansion rate is mostly large and posi¬ 
tive, leading to a bulk viscous pressure that reduces the 
effective pressure of the system and, consequently, slows 
down the acceleration of the fluid. 

As shown in Fig. 2, the calculations with bulk viscous 
pressure are indeed able to provide a good description of 
all the pT—integrated observables. The calculated average 
transverse momentum of pions, kaons, and protons are 
within the error bars of the ALICE measurements (SU] for 
most of the centrality classes considered. The pion and 
proton multiplicities measured by ALICE [50] are well de¬ 
scribed by the model, which however systematically over¬ 
predicts the multiplicity of kaons by ^ 10%. Finally, we 
see that the inclusion of bulk viscosity does not spoil the 
description of the flow harmonic coefficients U 2 , 3 , 4 { 2 } as 
a function of centrality. We note that the bulk viscosity 
reduces U 2 , 3 , 4 { 2 } by more than 10% but this effect is com¬ 
pensated by decreasing the shear viscosity over entropy 
density ratio from p/s = 0.16 to p/s = 0.095, leading to 
a very similar quality of description. Within this study, 
the inclusion of bulk viscosity can therefore reduce the 
value of shear viscosity extracted from data by almost 
50%. 

We now study py-differential observables within the 
best fit configuration including shear and bulk viscosi¬ 
ties. Figure 3 shows the pr-spectra of pions, kaons, and 
protons and U 2 , 3 , 4 { 2 }(pr) of charged hadrons for the 0- 
5% and 30-40% centrality classes. The solid lines cor¬ 
respond to the calculations with bulk and shear viscos¬ 
ity discussed above while the dashed lines correspond to 
the same calculations without the effect of hadronic re¬ 
scatterings. Note that the p^-spectra display reason¬ 
able agreement with the data which is in line with the 
good description of the multiplicity and (px) of pions. 
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FIG. 3: (Color online) Transverse momentum spectra (upper panels) of pions, kaons, and protons and harmonic flow coefficients 
(lower panels) as a function of the transverse momentum. Two centrality classes are considered: 0-5% (left panels) and 30- 
40% (right panels). The bands denote the statistical uncertainty of the calculation. The full and open symbols correspond to 
measurements by the ALICE [49] and CMS | 531I54| collaboration respectively, with bars denoting the experimental uncertainty. 


kaons, and protons, displayed in Fig. 2. The Vn{‘2}{pT) 
of charged hadrons shows more deviations from data, in 
particular the ALICE data [15], which is systematically 
smaller than the CMS measurement [53] [54] at high px- 
We find that hadronic re-scatterings have an almost 
negligible effect on pion spectra (the difference between 
the red dashed curve and the solid one is barely visible in 
the plot) and only affects the differential flow harmonics 
of charged hadrons at high pr- On the other hand, they 
play an important role in the description of kaon and, 
especially, proton spectra. Without taking into account 
all of these effects, it would not be possible to globally 
describe these observables. These findings are consistent 
with those from Refs. [sm^- 

4- Conclusions. In this letter, we discussed the effect 
of bulk viscous pressure on multiplicity, average trans¬ 
verse momentum, and azimuthal momentum anisotropy 
of charged hadrons using a state-of-the-art simulation 
of ultrarelativistic heavy ion collisions. It includes IP- 
Glasma initial conditions, which in combination with hy¬ 
drodynamics are known to provide a good description 
of the flow harmonic coefficients, and UrQMD, which 
models the hadronic re-scatterings that follow the fluid- 
dynamical evolution of the system. This fluid-dynamical 
evolution also considers several non-linear terms absent 
from several previous studies. The inclusion of bulk vis¬ 


cosity was found to have a large effect on the average 
transverse momentum of charged hadrons and on the el¬ 
liptic flow coefficient. In fact, when using the IP-Glasma 
initial conditions, the bulk viscosity is essential to de¬ 
scribe the pT-spectra of charged hadrons, and leads to 
a considerably better description of the data. A similar 
quality of description involving only shear viscosity could 
not be obtained in our current model. 

This work constitutes the first phenomenological in¬ 
vestigation which shows that the bulk viscosity of QGD 
matter is not small, at least around the phase transi¬ 
tion region. Our calculations suggest that C/s « 0.3 or 
larger around Tc- We also showed that the inclusion of 
bulk viscosity considerably modifies the optimum value 
of shear viscosity required to describe the data, reducing 
it by almost 50%. Therefore, the effects of bulk viscos¬ 
ity can not be neglected when extracting any transport 
coefficient from the data. The effects of bulk viscosity 
on ultracentral collisions, already briefly investigated in 
Ref. |46| . and on several other experimental observables 
will be the subject of future studies. 
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